home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / rng.exe / RAN.C < prev    next >
C/C++ Source or Header  |  1991-12-19  |  15KB  |  474 lines
  1. /* ran.c - functions to initialize random number generator and random number
  2. ** generator itself.
  3. **
  4. ** Module entry points
  5. ** -------------------
  6. **
  7. ** setup()  ----> N.B.  MUST be called before irand() or poiss()
  8. ** ran()    ----> macro defined in ran.h
  9. ** irand()
  10. ** poiss()
  11. ** binom()
  12. **
  13. **
  14. ** Either RAN_KNU or RAN_MAR must be defined before compilation, but
  15. ** only one can be defined.  If RAN_KNU is defined, a subtractive
  16. ** generator from volume two of Knuth is used.  If RAN_MAR is defined,
  17. ** a more complex generator discovered by Marsaglia is used.  Knuth's
  18. ** generator is about 50% faster than Marsaglia's, but its sequence is
  19. ** not as random.  For many purposes, Knuth's generator is good enough.
  20. ** It is much better than the generators normally supplied with C compilers,
  21. ** and it is often faster.
  22. **
  23. ** setup() must be called before irand() is called.  irand() does
  24. ** NOT check to see that the array it uses has been initialized.
  25. **
  26. ** irand() returns a uniformly distributed integer in [0..NO_NUMBERS)
  27. ** when NO_ZERO is not defined and in (0..NO_NUMBERS) when NO_ZERO is
  28. ** defined.
  29. **
  30. ** The macro ran() (defined in ran.h) returns a uniformly distributed
  31. ** random number in the range [0,1) when NO_ZERO is not defined and in
  32. ** (0,1) when NO_ZERO is defined.  The resolution is 1/CONV_REAL.
  33. **
  34. ** Strictly speaking poiss() should only be called when NO_ZERO is not
  35. ** defined.  But the error will be very small, since it will have an
  36. ** effect on the results in only 1 out of NO_NUMBERS calls, on average.
  37. **
  38. **
  39. ** References:
  40. **
  41. ** Knuth, D. E.  1981.  The Art of Computer Programming.  Vol. 2,
  42. **    Seminumerical Algorithms.  Addison-Wesley, Reading Mass.
  43. **
  44. **
  45. ** Author:  Kent E. Holsinger
  46. **
  47. ** Revision history
  48. ** ----------------
  49. **
  50. ** Version 1.0 -- 24 January 1990
  51. ** Original version of Knuth additive random number generator
  52. **
  53. ** Version 1.1 -- 24 January 1990
  54. ** ran() changed to macro operating on long
  55. **
  56. ** Version 2.1 -- 31 December 1990
  57. ** Added Marsaglia generator with conditional compilation
  58. **
  59. ** Version 2.3 -- 18 December 1991
  60. ** Converted both generators to use unsigned longs
  61. ** Editorial changes to internal comments
  62. **
  63. **
  64. */
  65.  
  66. #include <stdio.h>
  67. #include <math.h>
  68. #include "ran.h"
  69.  
  70.  
  71.                            /* constants used in random number generator */
  72. #ifdef RAN_KNU
  73. #define D      21
  74. #define LASTDIM      55
  75. #define FIRSTDIM  24
  76. #endif
  77. #ifdef RAN_MAR
  78. #define CD     7654321
  79. #define CM     16777213
  80. #define LASTDIM   98
  81. #define FIRSTDIM  34
  82. #endif
  83.  
  84.                            /*
  85.                            ** array of integers used in random number
  86.                            ** generator
  87.                            */
  88. typedef unsigned long INTARRAY[LASTDIM];
  89. typedef unsigned long *POINTER;
  90.  
  91.                            /* file used to store random number seed */
  92. FILE *randfile;
  93.  
  94.                            /* LOCAL VARIABLE DEFINITIONS */
  95. static INTARRAY x;         /* array used to generate random nos. */
  96. static POINTER j_ptr, k_ptr;  /* pointers into array */
  97.  
  98. #ifdef RAN_MAR
  99. static long c;
  100. #endif
  101. static char MOD_ID[] = "%W% %D%";
  102.  
  103.  
  104.                            /* LOCAL FUNCTION PROTOTYPES */
  105. #ifdef __STDC__
  106. static void rknuin(unsigned long m);
  107. static void rmarin(int ij, int kl);
  108. #else
  109. static void rknuin();
  110. static void rmarin();
  111. #endif
  112.  
  113.  
  114.  
  115. /***************************************************************************/
  116. /*                                                                         */
  117. /*                           MODULE ENTRY POINT                            */
  118. /*                                 setup()                                 */
  119. /*                                                                         */
  120. /***************************************************************************/
  121.                            /* Initialize the random number generator  */
  122. #ifdef __STDC__
  123. void setup(unsigned long *seed)
  124. #else
  125. void setup(seed)
  126.    unsigned long *seed;
  127. #endif
  128. {
  129. #ifdef RAN_MAR
  130.    unsigned long temp;     /* to write seed to file for next call */
  131.    int ij, kl;              /* seeds for Marsaglia generator */
  132. #endif
  133.  
  134.    randfile = fopen("RAND.INT", "r");
  135. #ifdef RAN_KNU
  136.    if (!randfile) {
  137.       printf("Enter a random number seed: ");
  138.       scanf("%lu", seed);
  139.    } else {
  140.       fscanf(randfile, "%lu", seed);
  141.       fclose(randfile);
  142.    }
  143.    do {
  144.       if (*seed == 0) {
  145.     fprintf(stderr, "Seed must be greater than zero\n");
  146.     fprintf(stderr, "Enter a random number seed: ");
  147.     scanf("%lu", seed);
  148.       }
  149.    } while (*seed == 0);
  150.    rknuin(*seed);
  151.    randfile = fopen("RAND.INT", "w");
  152.    fprintf(randfile, "%lu", x[LASTDIM - 1]);
  153.    fclose(randfile);
  154. #endif
  155. #ifdef RAN_MAR
  156.    if (!randfile) {
  157.       printf("Enter first random number seed:  ");
  158.       scanf("%d", &ij);
  159.       printf("Enter second random number seed: ");
  160.       scanf("%d", &kl);
  161.    } else {
  162.       fscanf(randfile, "%lu", seed);
  163.       fclose(randfile);
  164.       ij = (int) (*seed) % 31328;
  165.       kl = (int) ((*seed >> 16) % 30081);
  166.    }
  167.    if (ij < 0 || ij > 31328 || kl < 0 || kl > 30081) {
  168.       do {
  169.     fputs("The first random number seed must have a value between 0\
  170.  and 31328.\n", stderr);
  171.     fputs("The second seed must have a value between 0 and 30081.\n",
  172.           stderr);
  173.     printf("Enter first random number seed:  ");
  174.     scanf("%d", &ij);
  175.     printf("Enter second random number seed: ");
  176.     scanf("%d", &kl);
  177.       } while (ij < 0 || ij > 31328 || kl < 0 || kl > 30081);
  178.    }
  179.    rmarin(ij, kl);
  180.    temp = ((x[LASTDIM - 1] >> 16) % 30081) << 16;
  181.    temp += x[LASTDIM - 1] % 31328;
  182.    randfile = fopen("RAND.INT", "w");
  183.    fprintf(randfile, "%lu", temp);
  184.    fclose(randfile);
  185. #endif
  186. }
  187.  
  188.  
  189.  
  190. /***************************************************************************/
  191. /*                                                                         */
  192. /*                           MODULE ENTRY POINT                            */
  193. /*                                 irand()                                 */
  194. /*                                                                         */
  195. /***************************************************************************/
  196.                            /* Returns a random number in [0,NO_NUMBERS) */
  197.                            /* if NO_ZERO left undefined */
  198.                            /* Returns a random number in (0,NO_NUMBERS) */
  199.                            /* when NO_ZERO is defined */
  200. #ifdef RAN_KNU
  201.                            /* Knuth generator */
  202.  /*
  203.   * This is a direct translation of the subtractive generator (p. 171),
  204.   * which is based on the additive generator described on p. 27.  It has
  205.   * a period greater than 2^55 - 1.
  206.   *
  207.   * This implementation takes advantage of the fact that all elements of
  208.   * the array are themselves random numbers to return the element in the
  209.   * array currently pointed to.  Strict implementation of the algorithm
  210.   * would require definition of a temporary variable to hold the return
  211.   * value, which properly is *j_ptr before it is decremented
  212.   *
  213.   * On a Zenith Z-386 each call requires about 13 microseconds with
  214.   * NO_ZERO defined.  Without NO_ZERO defined each call requires about
  215.   * 11 microseconds.  (Results obtained using Microsoft C v6.00a,
  216.   * optimizing with /Ozx
  217.   */
  218. #ifdef __STDC__
  219. unsigned long irand(void)
  220. #else
  221. unsigned long irand()
  222. #endif
  223. {
  224.    *j_ptr -= *k_ptr;
  225.    j_ptr--;
  226.    k_ptr--;
  227.    if (j_ptr < x)
  228.       j_ptr = &x[LASTDIM - 1];
  229.    else if (k_ptr < x)
  230.       k_ptr = &x[LASTDIM - 1];
  231. #ifndef NO_ZERO
  232.    return *j_ptr;    /* simply return *j_ptr for U in [0,1) */
  233. #else
  234.    if (*j_ptr == 0)
  235.       return irand();
  236.    else
  237.       return *j_ptr;
  238. #endif
  239. }
  240. #endif
  241.  
  242. #ifdef RAN_MAR
  243.  /* Marsaglia generator */
  244.  /*
  245.   * C This random number generator originally appeared in "Toward a Universal
  246.   * C Random Number Generator" by George Marsaglia and Arif Zaman. C Florida
  247.   * State University Report: FSU-SCRI-87-50 (1987)
  248.   * C
  249.   * C It was later modified by F. James and pu